Harmonic oscillator eigenfunction expansions, quantum dots, and effective interactions 
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We give a thorough analysis of the convergence properties of the configuration-interaction method 
as applied to parabolic quantum dots among other systems, including a priori error estimates. The 
method converges slowly in general, and in order to overcome this, we propose to use an effective two- 
body interaction well-known from nuclear physics. Through numerical experiments we demonstrate 
a significant increase in accuracy of the configuration interaction method. 



I. INTRODUCTION 

The last two decades, an ever-increasing amount of 
research have been dedicated to understanding the elec- 
tronic structure of so-called quantum dotsi: semiconduc- 
tor structures confining from a few to several thousands 
electrons in spatial regions on the nanometre scale. In 
such calculations, one typically seeks a few of the lowest 
eigenenergies E k of the system Hamiltonian H and their 
corresponding eigenvectors iph, i.e., 



H%p k = E k ip k , 



(1.1) 



One of the most popular methods is the (full) configura- 
tion interaction method (CI) , where the many-body wave 
function is expanded in a basis of eigenfunctions of the 
harmonic oscillator (HO), and then necessarily truncated 
to give an approximation. In fact, the so-called curse of 
dimensionality implies that the number of degrees of free- 
dom available per particle is severely limited. It is clear, 
that an understanding of the properties of such basis ex- 
pansions is very important, as it is necessary for a priori 
error estimates of the calculations. Unfortunately, this is 
a neglected topic in the physics literature. 

In this article, we give a thorough analysis of the (full) 
configuration interaction (CI) method using HO expan- 
sions applied to parabolic quantum dots, and give practi- 
cal convergence estimates. It generalizes and refines the 
findings of a recent study of one-dimensional systems^ 
and is applicable to for example nuclear systems 3 - and 
quantum chemical calculations* as well. We demonstrate 
the estimates with calculations in the d = 2 dimensional 
case for N < 5 electrons, paralleling computations in the 
literature£&Z&9»i2. 

The main results are however somewhat discouraging. 
The expansion coefficients of typical eigenfunctions are 
shown to decay very slowly, limiting the accuracy of any 
practical method using HO basis functions. We therefore 
propose to use an effective two-body interaction to over- 
come, at least partially, the slow convergence rate. This 
is routinely used in nuclear physic s 3 ' 11 where the inter- 
particle forces are of a completely different, and basically 
unknown, nature. For electronic systems, however, the 
interaction is well-known and simpler to analyze, but ef- 
fective interactions of the present kind have not been ap- 
plied, at least to the author's knowledge. The modified 
method is seen to have convergence rates of at least one 



order of magnitude higher than the original CI method. 
An important point here is that the complexity of the 
CI calculations is not altered, as no extra non-zero ma- 
trix elements are introduced. All one needs is a relatively 
simple one-time calculation to produce the effective in- 
teraction matrix elements. 

The HO eigenfunctions are popular for several reasons. 
Many quantum systems, such as the quantum dot model 
considered here, are perturbed harmonic oscillators per 
se, so that the true eigenstates should be perturbations 
of the HO states. Moreover, the HO has many beautiful 
properties, such as complete separability of the Hamil- 
tonian, invariance under orthogonal coordinate changes, 
and thus easily computed eigenfunctions, so that comput- 
ing matrix elements of relevant operators becomes rela- 
tively simple. The HO eigenfunctions are defined on the 
whole of K d in which the particles live, so that trunca- 
tion of the domain is unnecessary. Indeed, this is one of 
the main problems with methods such as finite difference 
or finite element methods*^ On the other hand, the HO 
eigenfunctions are the only basis functions with all these 
properties. 

The article is organized as follows. In Sec.|TT]we discuss 
the harmonic oscillator and the the parabolic quantum 
dot model, including exact solutions for the N = 2 case. 
In Sec. IIII1 we give results for the approximation proper- 
ties of the Hermite functions in n dimensions, and thus 
also of many-body HO eigenfunctions. By approximation 
properties, we mean estimates on the error — Pip\\, 
where tp is any wave function and P projects onto a fi- 
nite subspace of HO eigenfunctions, i.e., the model space. 
Here, Pip is in fact the best approximation in the norm. 
The estimates will depend on analytic properties of tp, 
i.e., whether it is differentiable, and whether it falls of 
sufficiently fast at infinity. To our knowledge, these re- 
sults are not previously published. 

In Sec. lIVl we discuss the full configuration interaction 
method, using the results obtained in Sec. IIIII to obtain 
convergence estimates of the method as function of the 
model space size. We also briefly discuss the effective 
interaction utilized in the numerical calculations, which 
are presented in Sec. El We conclude with a discussion of 
the results, its consequences, and an outlook on further 
directions of research in Sec. [VTl 

We have also included an appendix with proofs of the 
formal propositions in Sec. IIIII 
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II. THE HARMONIC OSCILLATOR AND 
PARABOLIC QUANTUM DOTS 

A. The Harmonic Oscillator 

A spinless particle of mass m in an isotropic harmonic 
potential has Hamiltonian 

tfHO = -|V + imc 2 ||r1| 2 , (2.1) 

where f £ M. d is the particle's coordinates. By c hoosing 
proper energy and length units, i.e., hu and yjh/moj, 
respectively, the Hamiltonian becomes 



H W = -^V z + -\\f\\ 2 . 



(2.2) 



Huo can be written as a sum over d one-dimensional 
harmonic oscillators, viz, 



k=l 



k , 



(2.3) 



so that a complete specification of the HO eigenfunctions 
is given by 

$ ai , a2 ,...,a d (r) = ai (rO0a 2 (r 2 ) • • • (t>a d ( r d), ( 2 -4) 

where 4> ai (x), on = 0,1,... are one-dimensional HO 
eigenfunctions, also called Hermite functions. These are 
defined by 

<f> n (x) = {2 n nW 1 ' 2 r 1 ' 2 H n {x)e~ x2 l 2 1 n = 0, 1, . . . , 

(2-5) 

where the Hermite polynomials H„(x) are given by 



2 d n 

H n (x) = (-l)V — e- 



(2.6) 



The Hermite polynomials also obey the recurrence for- 
mula 



H n+ i(x) = 2xH n (x) - 2nH n -i(x) 



(2.7) 



with Ho(x) = 1 and Hi(x) — 2x. The Hermite poly- 
nomial H n (x) has n zeroes, and the Gaussian factor in 
4> n (x) will eventually subvert the polynomial for large 
\x\. Thus, qualitatively, the Hermite functions can be 
described as localized oscillations with n nodes and a 
Gaussian "tail" as x approaches ±oo. One can easily 
compute the quantum mechanical variance 

C°° 1 
(Ax) 2 := / x 2 <fi n {x) 2 dx = n + -, (2.8) 

J —oc 2 

showing that, loosely speaking, the width of the oscilla- 
tory region increases as (ri + 1/2) 1 / 2 . 

The functions ad defined in Eqn. (|2.4p are called 

d-dimensional Hermite functions. In the sequel, we will 
define a = (oil, ■ • • , ay) £ Xd for a tuple of non- negative 



integers, also called a multi-index; see Appendix I A II Us- 
ing multi-indices, we may write 

-1/2 



= 2Ha!^/ 2 H ai (n) ■ ■ ■ H ad {r d )e 



-||r|| 3 /2, 
(2.9) 

The eigenvalue of <fi n (x) is n + 1/2, so that the eigen- 
value of $ a (f) is 



e Q = 2 + N 



(2.10) 



i.e., a zero-point energy d/2 plus a non-negative inte- 
ger. We denote by \a\ the shell number of <i> Q , and 
the eigenspace 5^(18^) corresponding to the eigenvalue 
d/2 + r a shell. We define the shell-truncated Hilbert 
space V R {R d ) C L 2 (R d ) as 

R 

V R (R d ) := span{$ Q (f) | |a| < R) = 05 r (M d ), 

i.e., the subspace spanned by all Hermite functions with 
shell number less than or equal to R, or, equivalently, 
the direct sum of the shells up to and including R. The 
./V-body generalization of this space, to be discussed in 
Section fill Bl is a very common model space used in CI 
calculations. 

Since the Hermite functions constitute an orthonormal 
basis for L 2 (R d ), VB,(R d ) -> L 2 (R d ), in the sense that for 
every tp £ L 2 (R d ), lim^oo \\tp - Pip\\ = 0, where P is 
the orthogonal projector on Vn(M. d ). Strictly speaking, 
we should use a symbol like Pr or even Pfj(lR d ) for the 
projector. However, R and d will always be clear from the 
context, so we are deliberately sloppy to obtain a concise 
formulation. For the same reason, we will sometimes 
simply write V or Vr for the space VniW 1 ). 

An important fact is that since Huo is invariant under 
orthogonal spatial transformations (i.e., such transforma- 
tion conserve energy), so is each individual shell space. 
Hence, each shell S r (W i ), and also T']i(M. d ), is indepen- 
dent of the spatial coordinates chosen. 

For the case d = 1 each shell r is spanned by a single 
eigenfunction, namely <j) r (x). For d = 2, each shell r has 
degeneracy r + 1 , with eigenfunctions 

$ {s ,r- s ){7) = Mn)^s(r 2 ), 0<s<r. (2.12) 

The usual HO eigenfunctions used to construct many- 
body wave functions are not the Hermite functions 
^ai, -- .a d i however, but rather those obtained by utilizing 
the spherical symmetry of the HO. This gives a many- 
body basis diagonal in angular momentum. For d = 2 
we obtain the so-called Fock-Darwin orbitals given by 



2nl 



(n+\m\)\ 



1/2 gim6 



4:V) e - r2/2 . 



Here, n > is the nodal quantum number, counting the 
nodes of the radial part, and m is the azimuthal quantum 



3 



number. The eigenvalues are 



= 2n 



1. 



(2.14) 



Thus, R = 2n + \m\ is the shell number. By construc- 
tion, the Fock-Darwin orbitals are eigenfunctions of the 
angular momentum operator L z = —id/dg with eigen- 
value m. Of course, we may write as a linear 
combination of the Hermite functions $( s ,ii- s ), where 
< s < R — 2n + \m\, and vice versa. The actual 
choice of form of eigenfunctions is immaterial, as long as 
we may identify those belonging to a given shell. 

The space Vr = 4(M. 2 ) is illustrated in Fig. [T] using both 
Hermite functions and Fock-Darwin orbitals. 



B. Parabolic quantum dots 

We consider TV electrons confined in a harmonic oscil- 
lator in d dimensions. This is a very common model for a 
quantum dot. We comment, that modelling the quantum 
dot geometry by a perturbed harmonic oscillator is justi- 
fied by self-consistent calculations ; 13 ! 14 ' 15 and is a widely 
adopted assumptions^ ' 16 ' 17 ' 18 

The Hamiltonian of the quantum dot is given by 

H:=T + U, (2.15) 
where T is the many-body HO Hamiltonian, given by 



N 



(2.16) 



fc=i 



and U is the inter-electron Coulomb interactions. In di- 
mensionless units the interaction is given by, 



N 



N 



A 



(2.17) 



i<j 



The N electrons have coordinates f%, and the parame- 
ter A measures the strength of the interaction over the 
confinement of the HO, viz, 



Hlo \47T£oe / 



(2.18) 



where we recall that \Jhjmuj is the length unit. Typical 
values for GaAs semiconductors are close to A = 2, see 
for example Ref. [l8|. Increasing the trap size leads to 
a larger A, and the quantum dot then approaches the 
classical regime^ 



C. Exact solution for two electrons 

Before we discuss the approximation properties of the 
Hermite functions, it is instructive to consider the very 
simplest example of a two-electron parabolic quantum 



dot and the properties of the eigenfunctions, since this 
case admits analytical solutions for special values of A and 
is otherwise well understoo d 19 ' 20 ' 21 . Here, we consider 
d = 2 dimensions only, but the d = 3 case is similar. We 
note, that for N = 2 it is enough to study the spatial 
wave function, since it must be either symmetric (for the 
singlet 5 = spin state) or anti-symmetric (for the triplet 
5 = 1 spin states). The Hamiltonian (|2.15|) becomes 



H = -i(vf 



1 



A 



V|) + -(r 2 +r 2 ) + — 



ri2 



(2.19) 



where r 12 = \\rx — r^H an d T* = 



Introduce a set 



of scaled centre of mass coordinates given by R — (fx + 
7*2 VV^ and r = (fx — r%)I^Pl. This coordinate change 
is orthogonal and symmetric in R 4 . This leads to the 
separable Hamiltonian 



H 



1 



1 



(V^+V 2 fl ) + -(||r1 2 + || J R|| 2 ) + 



A 



2 , r n , 2 

= H W {R) + H, cl (f) 



V2\\r\\ 



A complete set of eigenfunctions of H can now be written 
on product form, viz, 



*(i?,r) = <P„ l ,( J R^(r) 



(2.20) 



The relative coordinate wave function tp{r) is an eigen- 
function of the relative coordinate Hamiltonian given by 



1 9 1 9 

-2 V ^:f 



A 

V2V 



(2.21) 



where r = ||r||. This Hamiltonian can be further sepa- 
rated using polar coordinates, yielding eigenfunctions on 
the form 



Am6 



(2.22) 



where |m| > is an integer and u n ^ m is an eigenfunction 
of the radial Hamiltonian given by 



1 d d \m\ 2 
-r- 



1 



A 



2rdr dr 2r 2 2 y/2r 



(2.23) 



By convention, n counts the nodes away from r = of 
u n ,m{ r )- Moreover, odd (even) m gives anti-symmetric 
(symmetric) wave functions \&(fi,7^). For any given \m\, 
it is quite easy to deduce that the special value A = 
•y/2|m| + 1 yields the eigenfunction 



-£>r |m| (a + r) e - r /2 , 



(2.24) 



where D and a are constants. The corresponding eigen- 
value of H r is E r = \m\ + 2, and E = 2n' + \m'\ + 1 + E r . 
Thus, the ground state (having m = m! = 0, n = n' = 0) 
for A = 1 is given by 

* (i?,r) = D{r + a)e- {r2+R2)/2 
D 



V2 



(r 12 + ^a) e -«+^/ 2 , 



4 



Ti — O n — 1 Ti — 2 n — 1 Ti — Q 

711 — —4 m — —2 m — 717 — 2 711 — 4 

71 — 71—1 71 — 1 71 — 

771 — —3 771 — —1 771 — 1 771 — 3 



a - (0, 4) q — (1, 3) ft — (2, 2) ft — (3, 1) ft - (4, 0) 



ft = (0,3) ft -(1,2) ft -(2,1) ft -(3,0) 



} shell <S 3 



n — n — l n — 

771 — — 2 771 — 771 — 2 



Unitarily equivalent 



a = (0,2) a = (1,1) a =(2,0) 



71 = 71 = 

771 = —1 771 = 1 

71 = 
771 = 



a = (0,1) a = (1,0) 



= (0,0) 



FIG. 1: Illustration of Pr=4,(R 2 ): (Left) Fock-Darwin orbitals. (Right) Hermite functions. Basis functions with equal HO 
energy are shown at same line. 



with D being a (new) normalization constant. 

Observe that this function has a cusp at r = 0, i.e., 
at the origin x = y = (where we have introduced 
Cartesian coordinates f = (x, y) for the relative coordi- 
nate). Indeed, the partial derivatives d x ^Po,o and dyipo.o 
are not continuous there, and 'J'o has no partial deriva- 
tives (in the distributional sense, see Appendix IA 2[) of 
second order. The cusp stems from the famous "cusp 
condition" which in simple terms states that, for a non- 
vanishing wave function at r\% = 0, the Coulomb diver- 
gence must be compensated by a similar divergence in the 
Laplaciani 22 i 23 This is only possible if the wave function 
has a cusp. 

On the other hand, the non-smooth function ^o(R,r) 
is to be expanded in the HO eigenfunctions, e.g., Fock- 
Darwin orbitals. (Recall, that the particular represen- 
tation for the HO eigenfunctions are immaterial - also 
whether we use lab coordinates 7*1,2 or centre-of-mass co- 
ordinates R and r, since the coordinate change is orthog- 
onal.) For to = 0, we have 

Or) = \[l L n (r*)e- ra ", (2.25) 
using the fact that these are independent of 9. Thus, 

oo 

* (r) = <Ofl)«o,o(r) = ^(R) E c «Or), 

(2.26) 

The functions <£>™ (r) are very smooth, as is seen by not- 
ing that L n (r ) = L n (x 2 + y 2 ) is a polynomial in x and 
y, while «o,o(t*) = u 0:0 {\/x 2 + y 2 ), so Eqn. (I2.26|) is ba- 
sically approximating a square root with a polynomial. 

Consider then a truncated expansion ^q^r € Pr(^- 2 ), 
such as the one obtained with the CI or coupled clus- 
ter method^ In general, this is different from Pr^q, 
which is the best approximation of the wave function in 
Vr(M ). In any case, this expansion, consisting of the 
R + 1 terms like those of Eqn. (|2.26[) is a very smooth 
function. Therefore, the cusp at r = cannot be well 
approximated. 



In Section MI C[ we will show that the smoothness 
properties of the wave function \P is equivalent to a cer- 
tain decay rate of the coefficients c n in Eqn. (|2.26[) as 
n — > oo. In this case, we will show that 

oo 

5> fc |c„| 2 < +00, (2.27) 

n=Q 

so that 

\c n \ =o(n-( k+1+ ^/ 2 ). (2.28) 

Here, k is the number of times VP may be differentiated 
weakly, i.e., "F e H k (M. 2 ), and e S [0,1) is a constant. 
For the function \l/o we have k = 1. This kind of esti- 
mate directly tells us that an approximation using only 
a few HO eigenfunctions necessarily will give an error 
depending directly on the smoothness k. 

We comment, that for higher |m| the eigenstates will 
still have cusps, albeit in the higher derivatives.— Indeed, 
we have weak derivatives of order |m| + 1, as can easily be 
deduced by operating on ipo,m with d x and d y . Moreover, 
recall that |m| = 1 is the 5=1 ground state, which 
then will have coefficients decaying faster than the 5 = 
ground state. Moreover, there will be excited states, i.e., 
states with |m| > 1, that also have more quickly decaying 
coefficients |c„|. This will be demonstrated numerically 
in Sec. El 

In fact, Hoffmann-Ostenhof et al£2- have shown that 
near ri2 = 0, for arbitrary A any local solution ^ of 
(H — E)ty = has the form 

*(0 = U\\ m P (1 + allfll) + 0(lieil" l+1 ), (2-29) 

where £ = (fi.f^) E M. 4 , and where P, deg(P) = to, 
is a hyper-spherical harmonic (on 5 3 ), and where a 
is a constant. This also generalizes to arbitrary N, 
cf. Sec. MIDI From this representation, it is manifest, 
that 'F S -ff m+1 (K 4 ), i.e., * has weak derivatives of or- 
der m+1. We discuss these results further in Sec. MI D"l 
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III. APPROXIMATION PROPERTIES OF 
HERMITE SERIES 

A. Hermite functions in one dimension 

In this section, we consider some formal mathemati- 
cal propositions whose proofs are given in Appendix I A 3( 
and discuss their importance for expansions in HO basis 
functions. 

The first proposition considers the one-dimensional 
case, and the second considers general, multidimensional 
expansions. The latter result has to the author's knowl- 
edge not been published previously. The treatment for 
one-dimensional Hermite functions is similar, but not 
equivalent to, that given by Boyd 2 ™ and Hille. 26 

We stress that the results are valid for any given wave- 
function - not only eigenfunctions of quantum dot Hamil- 
tonians - assuming only that the wavefunction decays ex- 
ponentially as |x| — > oo. In Appendix IA3| more general 
conditions are also considered. 

The results are stated in terms of weak differentiability 
of the wavefunction, which is a generalization of the clas- 
sical notion of a derivative. The space H k (M.) C L 2 (R) 
is roughly defined as the (square integrable) functions 
ip(x) having k (square integrable) derivatives d^ipix), 
< m < k. Correspondingly, the space H k (R n ) c 
L 2 (R") consists of the functions whose partial deriva- 
tives of total order < k are square integrable. For 
wavefunctions of electronic systems, it turns out that k 
times continuous differentiability implies k+1 times weak 
differentiability--. The order k of differentiability is not 
always known, but an upper or lower bound can often 
be found through analysis. It is however important, that 
the Coulomb singularity implies that k is finite. 

For the one-dimensional case, we have the following 
proposition: 

Proposition 1 (Approximation in one dimension) 

Let k > be a given integer. Let ip G L 2 (R) be expo- 
nentially decaying as \x\ — > oo and given by 

oo 

ip(x) = ^c n cf) n (x), (3.1) 

71=0 

where (p n {x) is given by Eqn. |2.5|) . Then ip 6 H k (R) if 
and only if 

oo 

$> fc |c„| 2 <oo. (3.2) 
n=0 

We notice that the latter implies that 

|C| =o(n-( fe+1 )/ 2 ), (3.3) 

which shows that the more ip(x) can be differentiated, the 
faster the coefficients will fall off as n — > oo. Moreover, 
let ip R = Pr^ = Sn=o c n4>n- Then 

/ oo \ 1/2 

u-m\ = [ E i c «i 2 > ( 3 - 4 ) 

\n=R+l ) 



which gives an estimate of how well a finite basis of Her- 
mite functions will approximate ip(x) in the norm. We 
already notice, that for low k = 2, which is typical, the 
coefficients fall off as o(n~ 3 / 2 ), which is rather slowly. 

In the general n-dimensional case, the wavefunction 
if) 6 L 2 (W n ) has an expansion in the n-dimensional Her- 
mite functions $ a (x), a € X„ given by 

a 

Qi ■ ■ -a n 

In order to obtain useful estimates on the error, we need 
to define the shell-weight p(R) by the overlap of ip( x ) 
with the single shell Sr, i.e., 

p(R) = \\P(SnW = E |c Q | 2 , (3.5) 

a, \cx\=R 

where P(Sr) is the projection onto the shell. Thus, 

oo 

||V|| 2 = X>0R). (3-6) 

R=0 

For the one-dimensional case, we of course have p(R) = 
\cr\ 2 - 

Proposition 2 (Approximation in n dimensions) 

Let %/} € i 2 (R") be exponentially decaying as \\x\\ — > oo 
and given by 

<P(x) =J2 c Mx). (3.7) 

a 

Then ip € H k (W l ) if and only if 

OO 

EH fc |c Q | 2 = E^W<+^- (3.8) 

a r—0 

Again, we notice that the latter implies that 

p(r) =o(r- (k+ V). (3.9) 

Moreover, for the shell-truncated Hilbert space Vr, the 
approximation error is given by 

/ oo XV 2 

ii(i-p)vii= E pW ■ ( 3J °) 

\r=R+l ) 

In applications, we often observe a decay of non- 
integral order, i.e., there exists an e G [0, 1) such that 
we observe 

p(r) =o(r-( k+1+ ^). (3.11) 

This does not, of course, contradict the results. To see 
this, we observe that if ip G H k (R n ) but ip H k+1 (R n ), 
then p(r) must decay at least as fast as o(r~( fc+1 )) but 
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not as fast as o(r~ k+2 ). Thus, the actual decay exponent 
can be anything inside the interval [k + 1, k + 2). 

Consider also the case where ip G H k (M. n ) for every 
k, i.e., we can differentiate it (weakly) as many times 
we like. Then p{r) decays faster than r~( k+1 \ for any 
k > 0, giving so-called exponential convergence of the 
Hermite series. Hence, functions that are best approxi- 
mated by Hermite series are rapidly decaying and very 
smooth functions ip. This would be the case for the quan- 
tum dot eigenfunctions if the inter-particle interactions 
were non-singular. 

B. Many-body wave functions 

We now discuss iV-body eigenfunctions of the HO in 
d dimensions, including spin, showing that we may iden- 
tify the expansion of a such with 2 N expansions in Her- 
mite functions in n = Nd dimensions, i.e., 2 N expan- 
sions in HO eigenfunctions of imagined spinless particles 
in n = Nd dimensions. Each expansion corresponds to a 
different spin configuration. 

Each particle k — 1, . . . , JV has both spatial degrees 
of freedom fk G M. d and a spin coordinate Tk G {±1}, 
corresponding to the z-projection S z = ±f of the elec- 
tron spin. The configuration space can thus be taken 
as two copies X of R d ; one for each spin value, i.e., 
X = R d x {±1} and Xk = {fk,Tk) € X are the coor- 
dinates of particle k. 

For a single particle with spin, the Hilbert space is now 
L 2 (X), with basis functions given by 

& i (x)=$ a (r) X *(T), (3.12) 

where i = i(a, a) is a new, generic index, and where Xa- 
is a basis function for the spinor space C 2 . 

Ignoring the Pauli principle for the moment, the N- 
body Hilbert space is now given by 

H(N) = L 2 (X) N = L 2 (R Nd ) $ (C 2 ) N , (3.13) 

i.e., each wavefunction ip £ Ti(N) is equivalent to 
2 N spin-component functions ip^ <= L 2 (R Nd ), a = 
(ax,-- - ,a N ) e {±1} N - We have 

i>{xi,-- ■ ,x N ) = ^2ip (a) (0x<r(T), £, = (n,- ■ ■ >f N ), 

(3-14) 

where r = (n, • • • , rjv), and where X<?( T ) = &o.t are basis 
functions for the iV-spinor space C 2 , being eigenfunc- 
tions for S z , i.e., corresponding to a given configuration 
of the N spins. 

The cr'th component function t/> (<t) (£) G L 2 (R Nd ) is 
a function of Nd variables, and by considering the Nd- 
dimensional HO as the sum of N HOs in d dimensions, 
it is easy to see that a basis for the L 2 (EL Nd ) is given by 
the functions 

^(O^a^i) (3.15) 



where £ = (fi, • ■ ■ , rjv), and where (3 = (a 1 , • • • , a N ) is 
an iVd-component multi-index. Correspondingly, a basis 
for the complete space H(N) = L 2 (X) N is given by the 
functions 

$ ll .., JV (C,T)E$ al (f 1 )...V(?Jv)x,W, (3.16) 

where ik — i(a , cr^). Notice, that the HO energy 
and hence the shell number only depends on (3 — 
(a\---,a N ). 

The functions ip^ may be expanded in the functions 
$0, i-e., 

^ {a Ho = £4^(0 

a 1 — a N 

and we define the cr'th shell weight p^ a \r) as before, i.e., 
pC">(r)= I^'l 2 - ( 3 - 17 ) 

|/3|=r 

We may then apply the analysis from Section IIII Al to 
each of the spin component functions, and note that the 
total shell-weight is 

since the shell number \/3\ = J2\ a k\ does not depend on 
the spin configuration of the basis function. 

Including the Pauli principle to accommodate proper 
wave-function symmetry does not change these consider- 
ations. The basis functions <&i 1 ...i N are anti-symmctrized 
to become Slater determinants Vp^...^ (see for example 
Ref. [27] for details), which is equivalent to consider the 
projection Has(^) — PasH(N) of the unsymmetrized 
space onto the antisymmetric subspace. Moreover, the 
projections Pr and Pas commute, so that the shell- 
truncated space is given by 

Vas,r = span I 4n : i x = i < ■ ■ ■ < i N , £ \a k \ < R 

{ k 

(3-19) 

which is precisely the computational basis used in many 
CI calculations. (See however also the discussion in Sec- 
tion[lV]) We stress, that, Vasm is independent of the ac- 
tual one-body HO eigenfunctions used. The shell-weight 
of -0 G Has(X) is now given by 

P(r) = £ (^ii-iN^) s \f3(ii-m)\,r, (3.20) 

ii—iN 

and 

R 

\\P R i>f =J2p(r). (3.21) 

r=Q 
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As should be clear now, studying approximation of 
Hermite functions in arbitrary dimensions automatically 
gives the corresponding many-body HO approximation 
properties, since the many-body cigcnfunctions can be 
seen as 2 N component functions, and since the shell- 
truncated Hilbert space transfers to a many-body setting 
in a natural way. 

C. Two electrons revisited 

We return to the exact solutions of the two-electron 
quantum dot considered in Sec. Ill CI Recall, that the 
wave functions were on the form 

fP(r,ff) = e" n6 f(r), (3.22) 

where f(r) decayed exponentially fast as r — > oo. Assume 
now, that tp 6 H k (R 2 ), i.e., that all partial derivatives of 
ip of order k exists in the weak sense, viz, 

^eL 2 (R 2 ), 0<j<k, (3.23) 

where x = rcos(9) and y = rsin(#). Then, by Lemma 0] 
in the Appendix, (aj, ) J : (aJ,) fc_J > € L 2 (R 2 ) for < j < k 
as well. 

The function ip(r, 0) was expanded in Fock-Darwin or- 
bitals, viz, 

oo 

V<r,0) = ^c„< D m M). (3.24) 

n=0 

Recall, that the shell number N for 3?^, was given by 
N = 2n + \m\. Thus, the shell- weight p(N) is in this case 
simply 

p(N) = |c (A r_| m | )/2 | 2 , N>\m\, (3.25) 
and p(N) = otherwise. From Prop. [2l we have 

oo 

N k p(N) < +oo, (3.26) 

N=\m\ 

which yields 

\c n \ = o(n-( fe + 1+e )/ 2 ), 0<e<l, (3.27) 
as claimed in Sec. Ill CI 

D. Smoothness properties of many-electron wave 
functions 

Let us mention some results, mainly due to Hoffmann- 
Ostenhof et a/. ) 22 ' 28 concerning smoothness of many- 
electron wave functions. Strictly speaking, their results 
are valid only in d = 3 spatial dimensions, since the 
Coulomb interaction in d — 2 dimensions fails to be a 
Kato potential, the definition of which is quite subtle and 



out of the scope for this article 28 . On the other hand, it is 
reasonable to assume that the results will still hold true, 
since the analytical results of the N = 2 case is very sim- 
ilar in the d — 2 and d — 3 cases: The eigenf unctions 
decay exponentially with the same cusp singularities at 
the originj 19 ! 20 

Consider the Schrodinger equation [H — E)ip(£) = 0, 
where f = (&,... ,£ Nd ) = (r u --- ,r N ) £ R Nd , and 
where is only assumed to be a solution locally. (A 
proper solution is of course also a local solution.) Recall, 
that i/j has 2 N spin-components ip^. Define a coalesce 
point £cp as a point where at least two particles coin- 
cide, i.e., fj = fi, j I. Away from the set of such 
points, tfj^ (£ ) is real analytic, since the interaction is 
real analytic there. Near a £cPj the wave function has 
the form 

^ (CT) (C + Ccp) = r k P(^/r)(l + ar) + 0(r k+1 ), (3.28) 

where r = ||£||, P is a hyper-spherical harmonic (on the 
sphere S Nd ~ 1 ) of degree k = fc(£cp)> and where a is a 
constant. It is immediately clear, that ip^"' (£ ) is k + 1 
times weakly differentiable in a neighborhood of £cp- 
However, at iiT-electron coalesce points, i.e., at points 
£cp where K different electrons coincide, the integer k 
may differ. Using exponential decay of a proper eigen- 
function, we have € H min ^ +1 (R Nd ). Hoffmann- 
Ostenhof et al. also showed, that symmetry restrictions 
on the spin-components due to the Pauli principle in- 
duces an increasing degree k of the hyper-spherical har- 
monic P, generating even higher order of smoothness. A 
general feature, is that the smoothness increases with the 
number of particles. 

However, their results in this direction are not general 
enough to ascertain the minimum of the values for k for 
a given wave function, although we feel rather sure that 
such an analysis is possible. Suffice it to say, that the 
results are clearly visible in the numerical calculations in 
Sec. El 

Another interesting direction of research has been un- 
dertaken by Yserentant, 2 - who showed that there are 
some very high order mixed partial derivatives at coa- 
lesce points. It seems unclear, though, if this can be 
exploited to improve the CI calculations further. 

IV. THE CONFIGURATION INTERACTION 
METHOD 

A. Convergence analysis using HO eigenfunction 
basis 

The basic problem is to determine a few eigenvalues 
and eigenfunctions of the Hamiltonian H in Eqn. (12.151) . 
i.e., 

Hifj k = E k i; k , k = 1, • • • , k max . (4.1) 

The CI method consists of approximating eigenvalues of 
H with those obtained by projecting the problem onto a 
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finite-dimensional subspace Tit C Tt(N). As such, it is 
an example of the Ritz-Galerkin variational methodi 30 i 31 
We comment, that the convergence of the Ritz-Galerkin 
method is not simply a consequence of the completeness 
of the basis functions^ We will analyze the CI method 
when the model space is given by 

H h = PrHas{N)=V r (N) 

= span | ^i lt ...,i N : 5^|a*| < r\ , 

used in Refs. [Tolf3^ . for example, although other spaces 
also are common. (We drop the subscript "AS" from now 
on.) The space 

M R (N) := span liw : max \a k \ < r\ , (4.2) 

i.e., a cut in the single-particle shell numbers (or energy) 
instead of the global shell number (or energy) is also 
common^ For obvious reasons, Vr(N) is often referred 
to as an "energy cut space", while A4r(N) is referred to 
as a "direct product space". 

As in Sec. MI A[ Pr is the orthogonal projector onto 
the model space Vr{N). We also define Qr = 1 — Pr 
as the projector onto the excluded space Vr(N) . The 
discrete eigenvalue problem is then 

(PrHPr)^^ = E h>k ip htk , k = 1, • • • , fc max . (4.3) 

The CI method becomes, in principle, exact as R — > oo. 
Indeed, a widely-used name for the CI method is "exact 
diagonalization," being somewhat a misnomer as only a 
very limited number of degrees of freedom per particle is 
achievable. 

It is clear that 

Vr{N) c Mr{N) c V nr {N), (4.4) 

so that studying the convergence in terms of Vr(N) is 
sufficient. In our numerical experiments we therefore fo- 
cus on the energy cut model space. A comparison be- 
tween the convergence of the two spaces is, on the other 
hand, an interesting topic for future research. 

Using the results in Refs. 13011331 for non-degenerate 
eigenvalues for simplicity, we obtain an estimate for the 
error in the numerical eigenvalue Eh as 

E h - E < [1 + v(R)](l + KX) (ip, Q R T^) , (4.5) 

where K is a constant, and where v(R) — ► as R — > oo. 
Using T§p = (Nd/2 + \/3\)$p and Eqn. IpTS]) , we obtain 

(^,QrT^)= {—+Ap{r). (4.6) 

r=R+l ^ ' 

Assume now, that ^ (ct) € H k (R Nd ) for all cr, so that 
according to Proposition [2l we will have 

oo 

^V'p(r) < +oo (4.7) 

r=0 



implying that rp(r) = o(r fe ). We then obtain, for k > 1, 
(V, (1 - P R )m = o(R-^) + o(R- k ). (4.8) 

For k — 1 (which is the worst case), we merely obtain 
convergence, (ip, (1 — Pr)Ti/j) — > as R — > oo. We as- 
sume, that R is sufficiently large, so that the o(R~ k ) term 
can be neglected. 

Again, we may observe a slight deviation from the de- 
cay, and we expect to observe eigenvalue errors on the 
form 

E h -E~ (l + ^A)i?- (fe - 1+£) , (4.9) 

where < e < 1. 

As for the eigenvector error — tp\\ (recall that tph ^ 
PriP), we mention that 

Hh - VII < [1 + r,(R)] [(1 + K\)(1>, (1 - P R )m\ 1/2 , 

(4.10) 

where t](R) — > as R — > oo . 

B. Effective interaction scheme 

Effective interactions have a long tradition in nuclear 
physics, where the bare nuclear interaction is basically 
unknown and highly singular, and where it must be 
renormalized and fitted to experimental data* 3 - In quan- 
tum chemistry and atomic physics, the Coulomb inter- 
action is of course well-known so there is no intrinsic 
need to formulate an effective interaction. However, in 
lieu of the in general low order of convergence implied 
by Eqn. (|4.9[) . we believe that HO-based calculations like 
the CI method in general may benefit from the use of 
effective interactions. 

A complete account of the effective interaction scheme 
outlined here is out of scope for the present article, but we 
refer to Refs. [^ ITlll34ll35ll36l for details as well as numerical 
algorithms. 

Recall, that the interaction is given by 

N N . 

i<j i<j h 1 j " 

a sum of fundamental two-body interactions. For the 
N = 2 problem we have in principle the exact solution, 
since the Hamiltonian (|2.19p can be reduced to a one- 
dimensional radial equation, e.g., the eigenproblem of H r 
defined in Eqn. (|2.21[) . This equation may be solved to 
arbitrarily high precision using various methods, for ex- 
ample using a basis expansion in generalized half-range 
Hermite functions^ In nuclear physics, a common ap- 
proach is to take the best two-body CI calculations avail- 
able, where R = O(10 3 ), as "exact" for this purpose. 

We now define the effective Hamiltonian for N — 2 as a 
Hermitian operator H Q g defined only within Vr(N = 2) 
that gives K = dim[Pfl(iV = 2)] exact eigenvalues E k 
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of H, and K approximate eigenvectors ip fi,k- Of course, 
there are infinitely many choices for the K eigenpairs, but 
by treating U = X/ri2 as a perturbation, and "following" 
the unperturbed HO eigenpairs (A = 0) through increas- 
ing values of A, one makes the eigenvalues unique 
The approximate eigenvectors tp e s,k € Vr(N — 2) are 
chosen by minimizing the distance to the exact eigenvec- 
tors ipk € H(N — 2) while retaining orthonormality. 35 
This uniquely defines H e g for the two-body system. In 
terms of matrices, we have 

H eS =Udiag(E 1 ,--- ,E K )U\ (4.12) 

where X and Y are unitary matrices defined as follows. 
Let U be the KxK matrix whose fc'th column is the coef- 
ficients of PRipk- Then the singular value decomposition 
of U can be written 

U = XYY\ (4.13) 
where £ is diagonal. Then, 

U:=XY^. (4.14) 

The columns of U are the projections PRipk "straight- 
ened out" to an orthonormal set. Eqn. (|4.12p is simply 
the spectral decomposition of H c g. Although different in 
form than most implementations in the literature (e.g., 
Ref. [ill ), it is equivalent. 

The effective two-body interaction C e s(i,j) is now 
given by 

Ceff(l, 2) := - PrTPr, (4.15) 

which is defined only within Vr(N = 2). 

The iV-body effective Hamiltonian is defined by 

N 

H cS := PrTPr + Ceff (*, j), (4-16) 

i<j 

where Pr projects onto 'Pr(N), and thus H e g is de- 
fined only within Vr{N). The diagonalization of H B g(N) 
is equivalent to a perturbation technique where a cer- 
tain class of diagrams is summed to infinite order in the 
full problem^ In implementations, (|4.12[) and (|4.16[) are 
treated in COM coordinates, utilizing block diagonality 
of both H and H e g, see Ref. [36| for details. 

We comment that unlike the bare Coulomb interac- 
tion, the effective two-body interaction C e g corresponds 
to a non-local potential due to the "straightening out" of 
truncated eigenvectors. 

Rigorous mathematical treatment of the convergence 
properties of the effective interaction is, to the author's 
knowledge, not available. Effective interactions have, 
however, enjoyed great success in the nuclear physics 
community, and we strongly believe that we soon will 
see sufficient proof of the improved accuracy with this 
method. Indeed, in Sec. [V] we see clear evidence of the 
accuracy boost when using an effective interaction. 



V. NUMERICAL RESULTS 
A. Code description 

We now present numerical results using the full 
configuration-interaction method for N = 2-5 electrons 
in d = 2 dimensions. We will use both the "bare" Hamil- 
tonian H = T + U and the effective Hamiltonian (|4.1G|) . 

Since the Hamiltonian commutes with angular momen- 
tum L z , the latter taking on eigenvalues M S Z, the 
Hamiltonian matrix is block diagonal. (Recall, that the 
Fock-Darwin orbitals are eigenstates of L z with 

eigenvalue m, so each Slater determinant has eigenvalue 
M = X)fcli TO fc-) Moreover, the calculations are done in a 
basis of joint eigenfunctions for total electron spin S 2 and 
its projection S z , as opposed to the Slater determinant 
basis used for convergence analysis. Such basis functions 
are simply linear combinations of Slater determinants 
within the same shell, and further reduce the dimension- 
ality of the Hamiltonian matrix^ The eigenfunctions of 
H are thus labeled with the total spin S = 0, 1, • • • , y 
for even N and S — \ , § , • • • , y for odd N, as well 
as the total angular momentum M = 0, l,---. (— M 
produce the same eigenvalues as M, by symmetry.) We 
thus split 'Pr(N) (or A4r(N)) into invariant subspaces 
Vr(N, M, S) (Mr(N, M, S)) and perform computations 
solely within these. 

The calculations were carried out with a code similar 
to that described by Rontani et al. in Ref. |8|. Table 
U shows comparisons of the present code with that of 
Table IV of Ref. 8 for various parameters using the model 
space M r (N, M, S) . Table U also shows the case A = 1, 
N = 2, M = 0, and 5 = 0, whose exact lowest eigenvalue 
is Eq = 3, cf. Sec. Ill CI We note that there are some 
discrepancies between the results in the last digits of the 
results of Ref. d. The spaces A4r(N) were identical in 
the two approaches, i.e., the number of basis functions 
and the number of non-zero matrix elements produced 
are cross-checked and identical. 

We have checked that the code also reproduces the 
results of Rcfs. 

mom 

using the Vr(N,M,S) spaces. 
Our code is described in detail elsewhere^ where it is 
also demonstrated that it reproduces the eigenvalues of 
an analytically solvable iV-particle system^ to machine 
precision. 



B. Experiments 

For the remainder, we only use the energy cut spaces 
Vr(N, M, S). Figure[5]shows the development of the low- 
est eigenvalue E = E (N, M, S) for N = 4, M = 0, 1, 2 
and S = as function of the shell truncation parame- 
ter R, using both Hamiltonians H and H c g. Apparently, 
the effective interaction eigenvalues provide estimates for 
the ground state eigenvalues that are better than the 
bare interaction eigenvalues. This effect is attenuated 
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TABLE I: Comparison of current code and Ref. Figures from the latter have varying number of significant digits. We include 
more digits from our own computation for reference 











R = 5 




R = 6 




R= 7 




N 


A 


M 


25 


Current 


Ref. 8 


Current 


Ref. 8 


Current 


Ref. 8 


2 


1 








3.013626 




3.011020 




3.009236 






2 








3.733598 


3.7338 


3.731057 


3.7312 


3.729324 


3.7295 






1 


2 


4.143592 


4.1437 


4.142946 


4.1431 


4.142581 


4.1427 


3 


2 


1 


1 


8.175035 


8.1755 


8.169913 




8.166708 


8.1671 




4 


1 


1 


11.04480 


11.046 


11.04338 




11.04254 


11.043 









3 


11.05428 


11.055 


11.05325 




11.05262 


11.053 


4 


6 








23.68944 


23.691 


23.65559 




23.64832 


23.650 






2 


4 


23.86769 


23.870 


23.80796 




23.80373 


23.805 


5 


2 





5 


21.15093 


21.15 


21.13414 


21.13 


21.12992 


21.13 




4 





5 


29.43528 


29.44 


29.30898 


29.31 


29.30251 


29.30 



14.6 




13.7 - 



13 ' 6 6 8 10 12 14 16 18 20 22 
R 

FIG. 2: Eigenvalues for N = 4, 3 = 0, X = 2 as function 
of R for H (solid) and H c g (dashed). M = 0, 1 and 2 are 
represented by squares, circles and stars, resp. 



with higher N, due to the fact that the two-electron ef- 
fective Coulomb interaction does not take into account 
three- and many-body effects which become substantial 
for higher N. 

We take the i? e ff-eigenvalues as "exact" and graph the 
relative error in Eq (N, M, S) as function of R on a loga- 
rithmic scale in Fig. [3l in anticipation of the relation 

ln(E h -E)KiC + alnR, a = -(k-l + e). (5.1) 

The graphs show straight lines for large R, while for small 
R there is a transient region of non-straight lines. For 
N = 5, however, A = 2 is too large a value to reach 
the linear regime for the range of R available, so in this 
case we chose to plot the corresponding error for the very 
small value A = 0.2, showing clear straight lines in the 
error. The slopes are more or less independent of A, as 
observed in different calculations. 



In Fig. 2] we show the corresponding graphs when using 
the effective Hamiltonian H c g. We estimate the relative 
error as before, leading to artifacts for the largest values 
of R due to the fact that there is a finite error in the best 
estimates for the eigenvalues. However, in all cases there 
are clear, linear regions, in which we estimate the slope a. 
In all cases, the slope can be seen to decrease by at least 
Aa w — 1 compared to Fig. [31 indicating that the effective 
interaction indeed accelerates the CI convergence by at 
least an order of magnitude. We also observe, that the 
relative errors are improved by an order of magnitude 
or more for the lowest values of R shown, indicating the 
gain in accuracy when using small model spaces with the 
effective interaction. 

Notice, that for symmetry reasons only even (odd) 
R for even (odd) M yields increases in basis size 
dim[P r(N , M, S)], so only these values are included in 
the plots. 

To overcome the limitations of the two-body effective 
interaction for higher N, an effective three-body interac- 
tion could be considered, and is hotly debated in the nu- 
clear physics community. (In nuclear physics, there are 
also more complicated three-body effective forces that 
need to be included4i) However, this will lead to a huge 
increase in memory consumption due to extra nonzero 
matrix elements. At the moment, there are no methods 
available that can generate the exact three-body effective 
interaction with sufficient precision. 

We stress, that the relative error decreases very slowly 
in general. It is a common misconception, that if a num- 
ber of digits of Eo (N, M, S) is unchanged between R and 
R + 2, then these digits have converged. This is not the 
case, as is easily seen from Fig. [3J Take for instance 
N = 4, M = and 5 = 0, and A = 2. For R = 14 
and R = 16 we have E = 13.84491 and E = 13.84153, 
respectively, which would give a relative error estimate of 
2.4 x 10~ 4 , while the correct relative error is 1.3 x 10~ 3 . 

The slopes in Fig. [3] vary greatly, showing that the 
eigenf unctions indeed have varying global smoothness, 
as predicted in Sec. IhTdI For (N,M,S) = (5,3,5/2), 
for example, a w —4.2, indicating that if) £ H 5 (R 10 ). It 
seems, that higher S gives higher k, as a rule of thumb. 
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Relative error for N=2, \ = 2 




Relative error for N=3, 1=2 



e — M=0, S=1 , a = -1 .2772 
* — M=0, S=3, a = -2.1716 
a — M=2, S=1, a = -1.3093 
A— M=2, S=3, a = -2.541 7 



18 20 




18 20 



Relative error for N=4, X = 2 



« — M=0, S=0, a = -1 .4233 
-* — M=0, S=4, a = -2.8023 
■e — M=3, S=2, a = -1.5327 
A— M=3, S=4, a = -3.2109 




Relative error for N=5, X = 0.2 



10 ' 



10" 



10 



-e — M=0, S=1 , a = -1 .51 50 
-* — M=0, S=5, a = -3.6563 
-B — M=3, S=1, a = -1.8159 
-A— M=3, S=5, a = -4.21 1 7 




18 20 



FIG. 3: Plot of relative error using the bare interaction for various N, M and S. Clear o(R a ) dependence in all cases. 



Intuitively, this is because the Pauli principle forces the 
wave function to be zero at coalesce points, thereby gen- 
crating smoothness. 

VI. DISCUSSION AND CONCLUSION 

We have studied approximation properties of Her- 
mite functions and harmonic oscillator eigenfunctions. 
This in turn allowed for a detailed convergence analy- 
sis of numerical methods such as the CI method for the 
parabolic quantum dot. Our main conclusion, is, that 
for wave functions ip € H k (M. n ) falling off exponentially 
as 1 1 a; 1 1 — > oo, the shell- weight function p(r) decays as 
p(r) = o(r~ k ~ 1 ). Applying this to the convergence the- 
ory of the Ritz-Galerkin method, we obtained the esti- 
mate (I4.5[) for the error in the eigenvalues. A complete 
characterization of the upper bound on the differentia- 



bility k, i.e., in ifj 1 -^ € H k , as well as a study of the 
constant K in Eqn. (14. 9f> , would complete our knowledge 
of the convergence of the CI calculations. 

We also demonstrated numerically, that the use of 
a two-body effective interaction accelerates the conver- 
gence by at least an order of magnitude, which shows 
that such a method should be used whenever possible. 
On the other hand, a rigorous mathematical study of the 
method is yet to come. Moreover, we have not investi- 
gated to what extent the increase in convergence is in- 
dependent of the interaction strength A. This, together 
with a study of the accuracy of the eigenvectors, is an 
obvious candidate for further investigation. 

The theory and ideas presented in this article should 
in principle be universally applicable. In fact, Figs. 1-3 
of Ref. |Tl| clearly indicates this, where the eigenvalues 
of 3 He as function of model space size are graphed both 
for bare and effective interactions, showing some of the 
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FIG. 4: Plot of relative error using an effective interaction for various N, M and S. Clear o(R a ) dependence in all cases, but 
notice artifacts when R is large, due to errors in most correct eigenvalues. The R = 5 case does not contain enough data to 
compute the slopes with sufficient accuracy. 



features we have discussed. 

Other interesting future studies would be a direct com- 
parison of the direct product model space Mr(N) and 
our energy cut model space Vr(N). Both techniques 
are common, but may have different numerical charac- 
teristics. Indeed, dim[A^i? \(N)] grows much quicker than 
dim[7- , it(iV)], while we are uncertain of whether the in- 
creased basis size yields a corresponding increased accu- 
racy. 

We have focused on the parabolic quantum dot, firstly 
because it requires relatively small matrices to be stored, 
due to conservation of angular momentum, but also be- 
cause it is a widely studied model. Our analysis is, how- 
ever, general, and applicable to other systems as well, 
e.g., quantum dots trapped in double-wells, finite wells, 
and so on. Indeed, by adding a one-body potential V 



to the Hamiltonian H = T + U we may model other 
geometries, as well as adding external fields i 42 ' 43 ' 44 
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APPENDIX A: MATHEMATICAL DETAILS 
1. Multi-indices 

A very handy tool for compact and unified notation 
when the dimension n of the underlying measure space 
W 1 is a parameter, are multi- indices. The set T n of multi- 
indices are defined as n-tuples of non-negative indices, 
viz, a — (ckj., • • • , a n ), where at > 0. 

We define several useful operations on multi-indices as 
follows. Let u be a formal vector of n symbols. Moreover, 
let 4>(0 = <K£i,£2, • • • , £n) be a function. Then, define 



a 

a±/3 
u a 



! = 



ai -+- «2 H h ot n 

a\\a2^- • ■ ■ a n \ 

(ai ± Pi, ■ ■ ■ ,a n ± Pn) 



(A.l) 
(A.2) 
(A.3) 
(A.4) 



gen gon 

Ct2 



T m (A.5) 



In Eqn. (|A.3[) . the result may not be a multi- index when 
we subtract two indices, but this will not be an issue for 
us. Notice, that Eqn. (|A.5|) is a mixed partial derivative 
of order |a|. Moreover, we say that a < ft if and only if 
o>j < 0j for all j. We define a = (j similarly. We also 



()x 



define "basis indices" ej by (ej)ji 
that we will often use the notation d x = — 
to simplify notation. Thus, 

consistent with Eqn. (|A.4|) . 



We comment, 
and d k = gf^ 



(A.6) 



2. Weak derivatives and Sobolev spaces 

We present a quick summary of weak derivatives and 
related concepts needed. The material is elementary and 
superficial, but probably unfamiliar to many readers, so 
we include it here. Many terms will be left undefined; 
if needed, the reader may consult standard texts, e.g., 
Ref . HE 

The space L 2 (R") is defined as 



L 2 (R n ) = \ 4> 



l^)| 2 d"£ < 



-co 



where the Lebesgue integral is more general than the 
Riemann "limit-of-small-boxes" integral. It is important 
that we identify two functions ip and %p\ differing only at 
a set Z £ K" of measure zero. Examples of such sets are 
points if n > 1, curves if n > 2, and so on, and countable 
unions of such. For example, the rationals constitute a 
set of measure zero in R. Under this assumption, L 2 (M. n ) 
becomes a Hilbert space with the inner product 



(A. 



where the asterisk denotes complex conjugation. 

The classical derivative is too limited a concept for the 
abstract theory of partial differential equations, including 
the Schrodinger equation. Let Cq° be the set of infinitely 
differentiable functions which are non-zero only in a ball 
of finite radius. Of course, Cg° C L 2 {« n ). Let -0 G 
L 2 (M n ), and let a S T n be a multi-index. If there exists 
awe L 2 (W l ) such that, for all <fi S Cfi°, 



(d^iOMOd^ = (-1) |Q| / <X£M£)d"£, (A.9) 



then d a ip = v G L 2 is said to be a weak derivative, or 
distributional derivative, of ip. In this way, the weak 
derivative is defined in an average sense, using integration 
by parts. 

The weak derivative is unique (up to redefinition on a 
set of measure zero), obeys the product rule, chain rule, 
etc. 

It is easily seen, that if i() has a classical derivative v S 
L 2 (R n ) it coincides with the weak derivative. Moreover, 
if the classical derivative is defined almost everywhere 
(i.e., everywhere except for a set of measure zero), then 
■0 has a weak derivative. 

The Sobolev space H k (R n ) is defined as the subset of 
L 2 (R n ) given by 

H k (R n ) ee {V G L 2 : <9 a G L 2 , Va G l n , \a\ < k) . 

(A.10) 

The Sobolev space is also a Hilbert space with the inner 
product 



(01,0 2 )= <d>l,9 a 02>, 



(A.H) 



|a|<fe 



and this is the main reason why one obtains a unified 
theory of partial differential equations using such spaces. 

The space H (K n ) for n > 1 is a big space - there 
are some exceptionally ill-behaved functions there, for 
example there are functions in H k that are unbounded 
on arbitrary small regions but still differentiable. (Her- 
mite series for such functions would still converge faster 
than, e.g., for a function with a jump discontinuity!) For 
our purposes, it is enough to realize that the Sobolev 
spaces offer exactly the notion of derivative we need in 
our analysis of the Hermite function expansions. 



3. Proofs of propositions 

We will now prove the propositions given in Sec. IIII1 
and also discuss these results on mathematical terms. 

Recall that the Hermite functions 4> n G L 2 (M), where 
n G No is a non-negative integer, are defined by 

</> n (x) = (2 n n\V^r 1,2 H n {x)e- x2 l 2 , (A.12) 
where H n (x) is the usual Hermite polynomial. 
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A well-known method for finding the eigenhmctions of 
i?HO hi onc dimension involves writing 

H„ = a*a + ^, (A.13) 
where the ladder operator a is given by 

a= -j=(x + d x ), (A.14) 
with Hcrmitian adjoint given by 

a f = -^(x-d x ). (A.15) 

The name "ladder operator" comes from the important 
formulas 

acp n {x) = y/n<j>n-i{x) (A. 16) 

d(f> n (x) = Vn+l<j> n+1 (x), (A.17) 

valid for all n. This can easily be proved by using the 
recurrence relation (|2.7p . By repeatedly acting on </> 
with (V we generate every Hermite function, viz, 

0„(x)=n!- 1 / 2 (at)> o (x). (A.18) 

As Hermite functions constitute a complete, orthonor- 
mal sequence in L 2 (R), any ip G L 2 (R) can be written as 
a series in Hermite functions, viz, 

oc 

ip(x) = ^c n (p n {x), (A.19) 

n=0 

where the coefficients c„ are uniquely determined by c n — 

An interesting fact is that the Hermite functions are 
also eigenfunctions of the Fourier transform with eigen- 
values (—*)", as can easily be proved by induction by 
observing firstly that the Fourier transform of (po( x ) is 
<f>o(k) itself, and secondly that the Fourier transform of 
is — ia> (acting on the variable y). It follows from 
completeness of the Hermite functions, that the Fourier 
transform defines a unitary operator on L 2 (R). 

We now make a simple observation, namely that 

{a}f<j> n {x) = Puinf^^+uix), (A.20) 

where Pk(n) = (n + k)\/n\ > is a polynomial of degree 
k for n > 0. Moreover Pf.(n + 1) > Pk{n) and Pk+i(n) > 
Pk(n). 

We now prove the following lemma: 

Lemma 1 (Hermite series in one dimension) 

Let ip G L 2 (R). Then 

1. atye L 2 (R) if and only ifaip G L 2 (R) if and only 
^ Y^=o n \ c n\ 2 < +oo, where c„ = (<p n ,ip) 

2. aty G L 2 (R) if and only ifxip, d x tp G L 2 (R) 



3. (at) fc +V G L 2 (R) implies (a)) k ip G L 2 (R) 

4. (at) k ip G i 2 (R) if and only if 

OO 

^n fc |c„| 2 < +oo. (A.21) 

n=Q 

5. (ot) k ip G i 2 (R) if and only if x j d^~ j ip G L 2 (R) for 
< j < k 

Proof: We have 

OO 

\\aH\\ 2 = £(n + l)|c„| 2 = ||^ 2 || + IMII 2 , (A.22) 

from which statement[T]follows. Statement[5]follows from 
the definition of , and that a'lp G L 2 implies atp G L 2 
(since ||a^|| < ||a^||), which again implies xip,d x ip G 
L 2 . Statement [3] follows from the monotone behaviour 
of Pk(n) as function of k. Statement 0] then follows. By 
iterating statement [5] and using [d x , x] = 1 and statement 
|3l statement [5] follows. 

The significance of the condition a*ip G L 2 (R) is that 
the coefficients c n of ip must decay faster than for a 
completely arbitrary wave function in L 2 (R). Moreover, 
aJip G L 2 (R) is the same as requiring d x ip G L 2 (R), 
and xtp G L 2 (R). Lemma [T] also generalizes this fact 
for (a') k ip G L 2 (R), giving successfully quicker decay of 
the coefficients. In all cases, the decay is expressed in 
an average sense, through a growing weight function in a 
sum, as in Eqn. (|A.21[) . Since the terms in the sum must 
converge to zero, this implies a pointwise faster decay, as 
stated in Eqn. (|A.25[) below. 

We comment here, that the partial derivatives must be 
understood in the weak, or distributional, sense: Even 
though ip may not be everywhere differentiable in the 
ordinary sense, it may have a weak derivative. For ex- 
ample, if the classical derivative exists everywhere except 
for a countable set (and if it is in L 2 ), this is the weak 
derivative. Moreover, if this derivative has a jump dis- 
continuity, there are no higher order weak derivatives. 

Loosely speaking, since xip(x) G L 2 4=> d y ip{y) G L 2 , 
where tp(y) is the Fourier transform, point [2] of Lemma 
[T]is a combined smoothness condition on ip(x) and ip(y). 
Point [5] is a generalization to higher derivatives, but is 
difficult to check in general for an arbitrary ip. On the 
other hand, it is well-known that the eigenfunctions of 
many Hamiltonians of interest, such as the quantum dot 
Hamiltonian (|2.15p . decay exponentially fast as \x\ — ► oo. 
For such exponentially decaying functions over R , x k tp G 
i 2 (R) for all k > 0, i.e., ip(y) is infinitely differentiable. 
We then have the following lemma: 

Lemma 2 (Exponential decay in ID) 

Assume that x k ip G L 2 (R) for all k > 0. Then a sufficient 
criterion for (o)) m ip G L 2 (R) is d™ip G L 2 (R), i.e., ip G 
H m (R). In fact, x k d x n 'ip G L 2 for all m' < m and all 
k>0. 



Proof: We prove the Proposition inductively. We note 
that, since d™ip £ L 2 implies € L 2 , the proposi- 

tion holds for m — 1 if it holds for a given m. Moreover, 
it holds trivially for m = 1. 

Assume then, that it holds for a given m, i.e., that 
ip £ H m implies x k d x n ~~i?p £ L 2 for 1 < j < m and for 
all k (so that, in particular, (a^) m i/j £ L 2 ). It remains to 
prove, that ip £ H m+1 implies x k d x n ip £ L 2 , since then 
(a^) m+1 ?/> £ L 2 by statement [5]of Lemma[TJ We compute 
the norm and use integration by parts, viz, 

\\x k d^\\ 2 = J^x 2k d k x r(x)d k x ^{ X ) 

= -2fc(9™V',^ 2fc ~ 1 a" l ~V> 
~<9™ + V,^ 2fc 5™-V> < +oo. 

The boundary terms vanish. Therefore, x k d x n ip £ L 2 for 
all k, and the proof is complete. § 

The proposition states that for the subset of L 2 (R) 
consisting of exponentially decaying functions, the ap- 
proximation properties of the Hermite functions will only 
depend on the smoothness properties of ip. Moreover, the 
derivatives up to the penultimate order decay exponen- 
tially as well. (The highest order derivative may decay 
much slower.) 

From Lemmas [1] and [5] we extract the following impor- 
tant characterization of the approximating properties of 
Hermite functions in d = 1 dimensions: 

Proposition 1 (Approximation in one dimension) 

Let k > be a given integer. Let ip £ L 2 (R) be given by 

oo 

C-n 

). (A.23) 

n=0 

Then ip £ H k (R) if and only if 

OO 

^V|c„| 2 <(5o. (A.24) 

The latter implies that 

\c n \ =o(n-( fe+1 )/ 2 ). (A.25) 

Let ip R = P R ip = J2n=o c n4>n- Then 

/ oo \ 1/2 

M>"lM = E l c «! 2 ■ (A.26) 

\n=R+l ) 

This is the central result for Hermite scries approxima- 
tion in L 2 (R 1 ). Observe that Eqn. (|A.25j) implies that 
the error \\ip — ip R \\ can easily be estimated. See also 
Prop. [Hand comments thereafter. 

Now a word on pointwisc convergence of the Her- 
mite series. As the Hermite functions are uniformly 
bounded^ viz, 

\<f> n (x)\ < 0.816 VxeR, (A.27) 
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the pointwise error in tp R is bounded by 

oo 

- i> R (x)\ < 0.816 £ |c„|. (A.28) 

n=R+l 

Hence, if the sum on the right hand side is finite, the con- 
vergence is uniform. If the coefficients c n decay rapidly 
enough, both errors can be estimated by the dominating 
neglected coefficients. 

We now consider expansions of functions in L 2 (R"). 
To stress that M™ may be other than the configura- 
tion space of a single particle, we use the notation x — 
(xi, • • • , x n ) £ W 1 instead of r £ R d . For N electrons in 
d spatial dimensions, n = Nd. 

Recall that the Hermite functions over R™ are indexed 
by multi-indices a = {a\, • • • , a n ) £ X n , and that 

$ a (x) = 4> ai {xi) (x n ). (A.29) 

Now, to each spatial coordinate Xk define the ladder op- 
erators ctfc = (xk + 9fc)/v2. These obey [oj, a&] = and 
[dj, a k ] = <5j.fc, as can easily be verified. Let a be a formal 
vector of the ladder operators, viz, 

a = (ai,a 2 , . . . ,a n ). (A. 30) 

For the first Hermite function, we have 

$ (ar) =Tr~ n/4 e-^ 2 / 2 . (A.31) 

By using Eqn. (|A.18|) and Eqn. (|A.4|) . we may generate 
all other Hermite functions, viz, 

$ a (x) = a!- 1 / 2 ( a t)«$ ( x ). (A.32) 

Given two multi-indices a and (3, we define the poly- 
nomial P a (f3) by 

P {a) = ^±$l = f[p {aj): (A.33) 
i=i 

where P is defined for integers as before. 

Since the Hermite functions $ Q constitute a basis, any 
if) £ L 2 (W l ) can be expanded as 

a a 

where the sum is to be taken over all multi-indices a £ 
T n . Now, let (3 £ l n be arbitrary. By using Eqn. (|A. 1 7|) 
in each spatial direction we compute the action of (a^)^ 
on tp: 

(atyV = (a\r^---(al)^J2 c ^ 

a 

^ i r K + /3 fc )! 1/2 „ 

a fc=l k ' 
a 
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Similarly, by using Eqn. (|A.16[) we obtain 

a 

a fe=l K ' 

= ^Tc^P^a) 1 ^, (A.36) 

a 

Computing the square norm gives 

||(at)^|| 2 = ^P^(a)|c Q | 2 . (A.37) 

a 

and 

||a^|| 2 =£p^(a)|c a+/3 | 2 . (A.38) 

a 

The polynomial Pp(ot) > for all (3,a E I n , and Pp{a + 
a') > Ppfa) for all non-zero multi-indices a' ^ 0. There- 
fore, if (at)<V G L 2 (R") then a^ifj E L 2 {R n ). However, 
the converse is not true for n > 1 dimensions, as the norm 
in Eqn. ((A. 38(1 is independent of infinitely many coeffi- 
cients c Q , while Eqn. (|A.37|) is not. (This should be con- 
trasted with the one-dimensional case, where aip E L 2 (R) 
was equivalent to a'ip E L 2 (R).) On the other hand, as in 
the n = 1 case, the condition alip E L 2 (R") is equivalent 
to the conditions x k ip E L 2 (R n ) and d k ip E L 2 (R n ). 

We are in position to formulate a straightforward gen- 
eralization of Lemma [TJ The proof is easy, so we omit it. 

Lemma 3 (General Hermite expansions) 

Let ip E L 2 (R n ), with coefficients c a as in Eqn. &A.34\) , 
and let (3 E T n be arbitrary. Assume (a t ) /3 V G L 2 (R n ). 
Then (aJ)P'ip E L 2 (R n ) and aP'ip E L 2 (R") for all (3' < 
(3. Moreover, the following points are equivalent: 

1. (a^tp g L 2 (R") 

2. For all multi-indices 7 < j3, x^d^'^tp E L 2 (R n ). 
3- Logical 2 < +00 

We observe, that as we obtained for n = 1, condition 
[5] is a combined decay and smoothness condition on ip, 
and that this can be expressed as a decay-condition on 
the coefficients of ip in the Hermite basis by [3l 

Exponential decay of ip E L 2 (R n ) as ||aj|| — > 00 im- 
plies that that x^ip E L 2 (R n ) for all 7 6 I„. We now 
generalize Lemma [2] to the n-dimensional case. 

Lemma 4 (Exponentially decaying functions) 

Assume that ip E L 2 (R n ) is such that for all 7 6 I„, 
x^ip G L 2 (R"). Then, a sufficient criterion for (a^^tp G 
L 2 (R n ) is d^ip E L 2 (R n ). Moreover, for all fi < (3, we 
have x^dP-^ip E L 2 {R n ) for all 7 G l n such that 7fe = 
whenever fi k = 0, i.e., the partial derivatives of lower 
order than (3 decay exponentially in the directions where 
the differentiation order is lower. 



Proof: The proof is a straightforward application of 
the n — 1 case in an inductive proof, together with the 
following elementary fact concerning weak derivatives: If 
1 < j < k < n, and if Xjip(x) and dkip{x) are in L 2 (R"), 
then, by the product rule, dk(xjip(x)) — Xj(dkip(x)) E 
L 2 (R n ). Notice, that Lemma [2] trivially generalizes to a 
single index in n dimensions, i.e., to j3 = Pk&k> since the 
integration by parts formula used is valid in R n as well. 
Similarly, the present Lemma is valid in n — 1 dimensions 
if it holds in n dimensions, as it must be valid for (3 = 

(0,0a,- •• ,Pn). 

Assume that our statement holds for n — 1 dimensions. 
We must prove that it then holds in n dimensions. As- 
sume then, that d^ip E L 2 {R n ). Let = d^tp E L 2 . 
Moreover, d^<p E L 2 . Since ip is exponentially decaying, 
and by the product rule, x'^cp E L 2 for all 71 > 0. By 
Lemma xj'd^'^cp G L 2 for all 71 and < /Zi < ft. 
Thus, xl 1 d /3 - ei,11 i[> E L 2 . Thus, the result holds as long 
as /i = ei/ii; or equivalently p — e k \Xk for any k. To 
apply induction, let \ = xj 1 di 1 ~ ,>1 ip E L 2 . Note that 
d^x e L 2 and x^ X € L 2 for all 7 = (0,72, • • • ,7„). But 
by the induction hypothesis, x^d^~^x S L 2 for all p, < (3 
and all 7 such that 7& = if fif. = 0. This yields, using 
the product rule, that x 1 d^~^ip E L 2 for all fx < j3 and 
all 7 such that 7^ = if fi k = 0, which was the hypoth- 
esis for n dimensions, and the proof is complete. Notice, 
that we have proved that (a')"tp E L 2 as a by-product. 



In order to generate a simple and useful result for 
approximation in n dimensions, we consider the case 
where ip decays exponentially, and ip € H k (R n ), i.e., 
d^ip(x) E L 2 (R n ) for all (3 G 2 n w ith = k. In this 
case, we may also generalize Eqn. (|A.25[) . For this, we 
consider the shell-weight p(r) defined by 

p(r)= J2 l c «l 2 < (A.39) 

a £T n , I ot\ =r 

where c Q = ($ a ,ip). Then, ||-0|| 2 = J2^LoP( r )- More- 
over, if P projects onto the shell-truncated Hilbert space 
P fl (R n ), then 

||P^|| 2 = 5>(r). (A.40) 

r=0 

Proposition 2 (Approximation in n dimensions) 

Let ip E L 2 (R") be exponentially decaying and given by 

i)(x) =Y,cMx)- (A.41) 

a 

Then ip E H k (R n ), k>0, if and only if 

OO 

\<A k K? = r "p^ < +°°- ( A - 42 ) 

a r— 

The latter implies that 

p{r)=o{r- { - k+1 ^). (A.43) 
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Moreover, for the shell-truncated Hilhert space Vr, the 
approximation error is given by 

/ oo \ 1/2 

||(1-P)V||= E ■ (A.44) 

\r=R+l / 

Proof: The only non-trivial part of the proof concerns 
Eqn. (|A.42[) . Since tp is exponentially decaying and since 
ip e H k if and only if d^ip G L 2 for all (3, \[3\ < k, we 
know that ^2 a a^\c a \ 2 < Too for all f), \[3\ < k. Since 



\a\ k is a polynomial of order k with terms of type apa? , 
d/3 > and |/3| = fc, we have 

E|a| fe |c a | 2 = E a^E^I^I 2 < +°°- ( A - 45 ) 

a /3, \0\=k a 

On the other hand, since ap > and the sum over 
j3 has hnitely many terms, ^ Q |a| fc |c a | < +oo implies 
X^c^lcal 2 < +oo for all /?, |/3| = fc, and thus ip G H k 
since tp was exponentially decaying. 
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